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Correlations in the motion of reptating polymers in a melt are investigated by means of Monte 
f*^ ■ Carlo simulations of the three dimensional slithering snake version of the bond-fluctuation model. 

' Surprisingly, the slithering snake dynamics becomes inconsistent with classical reptation predictions 

at high chain overlap (created either by chain length N or by the volume fraction <f) of occupied 
lattice sites), where the relaxation times increase much faster than expected. This is due to the 
anomalous curvilinear diffusion in a finite time window whose upper bound t+(N) is set by the 
density of chain ends <j>/N. Density fluctuations created by passing chain ends allow a reference 
polymer to break out of the local cage of immobile obstacles created by neighboring chains. The 
' dynamics of dense solutions of "snakes" at t <C r+ is identical to that of a benchmark system where 

all chains but one are frozen. We demonstrate that the subdiffusive dynamical regime is caused by 
the slow creeping of a chain out of its correlation hole. Our results are in good qualitative agreement 
«4 ■ with the activated reptation scheme proposed recently by Semenov and Rubinstein [Eur. Phys. J. B, 

O ' 1 (1998) 87]. Additionally, we briefly comment on the relevance of local relaxation pathways within 

a slithering snake scheme. Our preliminary results suggest that a judicious choice of the ratio of 
local to slithering snake moves is crucial to equilibrate a melt of long chains efficiently. 
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PACS numbers: 61.25.Hq Macromolecular and polymer solutions; polymer melts; swelling — 83.10.Kn Rep- 
, tation and tube theories - 83.80.Sg Polymer melts 
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JL ■ I. INTRODUCTION 

Background. The molecular-level description of the dynamics of non-dilute polymer liquids remains a challenging 
problem of great fundamental interest as well as practical importance 0. In the last 30 years, attention has largely 
focused on the applicability of the reptation concept introduced by de Gennes and Edwards 0, • This approach 
postulates that individual chains, constrained by their neighbors, move primarily along their own contours in a 
^ ■ snakelike fashion. Reptation-based models have been very successful in reconciling a wide range of experimental 
observations and are broadly in agreement with more recent computational results as reviewed in 0, All 
evidence amounts to the existence of an additional and chain length independent length scale. This length scale is 
called the "tube diameter" d e ~ N£ with jV e (</>, Z p ) the associated "entanglement length", <f> the molecular density, l p 
the persistence length and v the Flory exponent of the polymer chain. The density and persistence length dependence 
' of d e have been probed by various numerical methods [y, Q, S H3 , indicating that it is proportional to the excluded 
volume blob size £(</>, l p ). While, strictly speaking, not having the (more prestigious) status of a 'scientific theory', the 
"j***; i reptation concept is widely (though not generally) seen as a physically appealing starting point any further modeling 
& \ approach has to compete with. 

First motivation. The simulation of entangled polymer melts is a challenging task and it is a matter of debate if 
the (suspected) reptating motion of polymer chains has already been reproduced numerically in a sufficiently large 
simulation box Without any doubt, it will be near to impossible to simulate within the near future by means 
of locally realistic dynamics industrially relevant problems, such as the rheology of polymer melts with fillers. It is, 
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therefore, tempting to accept the snakelike motion as a basic ingredient in a coarse-grained algorithm and to verify if 



a meaningful dynamic interpretation of the well-known "slithering snake algorithm" (SSA) |12|, [13J, [1J| can be given. 
This algorithm is sketched in Fig. ^ It has been extensively used for three decades in simulations of static properties 
because of its efficiency in exploring phase space. Clearly, it would be extremely useful if a (simple and transparent) 



mapping of this algorithm on its corresponding local dynamical scheme could be achieved. 
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Self- consistency of reptation. While most of the current computational work appears to focus on the demonstration 
of the curvilinear chain motion [HI lil ITcL fTTL Hfij. here we take this key reptation hypothesis for granted. Therefore, 
we are not interested in the Rouse-like dynamics at short times and for chain length smaller than N e . For the SSA 
we expect N c 1 so that it should be possible to explore the dynamics of the so-called 'primitive path' What 
we do question, however, is an additional (and arguably less profound) mean-field assumption made in the original 
versions of reptation which supposes the uncorrelated motion of neighboring "snakes". 1 In the following we refer to 
this assumption as the "uncorrelated reptation hypothesis" (URH). One of its immediate consequences is, for instance, 
the linear time dependence of the mean-square displacements (MSD) of the center of mass or of the curvilinear motion 
of an inner chain monomer for the snake dynamics. The verification of the influence of the interchain correlations on 
reptating snakes is the second motivation of this study. 

Internal modes, tube fluctuations. In recent years a number of theoretical models and numerical schemes jjj ITfiL 
ITtI Il8l| have been developed in order to explain the observed molecular weight dependence of the zero-shear viscosity 
r\ ~ 7\r3-4±o.2 £ Qr TV > iV e 0,0. The main idea of these theories is to account for the internal modes of polymer 
dynamics inside the reptation tube ( "tube-length fluctuations" ) or to introduce a finite lifetime of the tube ( "double 
reptation" 01) without addressing, however, the possibility of correlated chain motion. While these are undoubtedly 
valiant efforts, they are of no concern here, as we are, to stress it again, interested in the physics at much larger 
molecular weights and correspondingly longer time scales. 

Coupling between dynamics and density fluctuations. The possibility of a coupling between primitive path dynam- 
ics and the fluctuating excluded volume interactions of neighboring chains was first investigated by Deutsch He 
suggested that the reptation of a chain tail from its original tube to a new environment implies a finite free energy 
penalty F(s) proportional to the curvilinear length of the tail s. Moreover, it was assumed that this additional acti- 
vation energy of the snake can relax as soon as the spatial size of the escaped portion becomes of the order of the 
mean distance between chain ends d cn d ~ A 1 / 3 . Taking into account the Gaussian chain statistics in a dense polymer 
solution, Deutsch obtained an exponential increase r ~ 77 ~ N 3 exp(A 2 / 3 ) for the relaxation times and viscosities in 
concentrated polymer solutions with extremely large molecular masses. More recently, these ideas have been elabo- 
rated further by Semenov and Rubinstein [20l l2ll l22| and implemented in a slig htly different and thermodynamically 
fully consistent way. They predict that the URH breaks down for N N 3 22]. Since we expect N e to be of order 1 
for the SSA, it should be possible to test this "activated reptation!^ hypothesis (ARH) in our simulations. 

Brief summary of our work. We have harnessed a SSA of the bond-fluctuation lattice model and compared its 
dynamics for consistency with the URH of the classical reptation model. We find that the dynamics rapidly slows 
down with chain overlap and becomes strongly subdiffusive in a finite time interval t_ < t « t + whose bounds are 
determined by the blob size £ and the mean distance between chain ends d G nd- This slowing down is at variance with the 
URH and indicates strong correlations between neighboring chains. To corroborate our results for the SSA dynamics 
in an annealed molecular field, referred to as "annealed SSA" (aSSA) hereafter, we have considered a different version 
of the SSA where all chains but one have been quenched. We will refer to this version as the "quenched SSA" (qSSA). 
The comparison of both dynamical schemes demonstrates the mutual influence of the correlation hole retaining the 
snakes by means of an effective force and the fluctuating chain ends which allow the snakes to sneak out of the local 
cage of immobile monomer sites. 

Outline. This work is presented as follows. Our numerical scheme is briefly recapitulated in the next section. The 
parameter range and some technical points concerning the configurations and measurements are presented in Sec. lIIII 
The relevant analytical predictions according to the URH and ARH are summarized in Sec. II VI The numerical results 
of this paper are presented in Sec. IVl where we characterize and compare aSSA and qSSA dynamics. The implications 
and possible extensions of our work are discussed in the final Sec. I VII There, we speculate on the importance of local 
motion to obtain a fully consistent mapping of a slithering snake scheme on a locally realistic model. 



II. SIMULATION MODEL 



Bond-fluctuation model. In this investigation we use the extensively studied 'bond-fluctuation model' (BFM) 
[7L I23I l24j . a lattice Monte Carlo scheme introduced by Carmesin and Kremer [23 |. Figure Q depicts a two dimensional 
projection of both the original BFM with local monomer dynamics (upper chain) and the BFM with global snakelike 
moves (lower chain) . These versions have been discussed elsewhere in detail [2^, 0, HH and we concentrate here on 
some salient features. 



We will often use the term "snake" as a short form for "chains moving according to the slithering snake algorithm" . 



In the BFM each monomer occupies the eight corners of the unit cube on a simple cubic lattice. The bonds between 
adjacent monomers can vary in length and direction, subject only to excluded volume constraints and entanglement 
restrictions. Basically, due to the larger number of degrees of freedom — 87 angles between bonds are possible — the 
BFM does not suffer from the same ergodici ty p roblems that classical Monte Carlo schemes of self-avoiding random 
chains on simple cubic lattices are prone to (We come back to this point in Sec. IV El ) In the original BFM 

version only local jumps to nearest neighbor sites are allowed. The excluded volume then automatically forbids any 
intersection of bonds. 

Slithering snake schemes: aSSA and qSSA. If snakelike moves are considered, the uncrossability of bonds has to 
be checked explicitly to avoid configurations which cannot be attained or unraveled by the original local scheme [27| ■ 
The snake move along the existing backbone of a chain is equivalent to removing an end monomer and connecting 
it to the other end of the chain, leaving the lattice positions of the middle monomers unchanged. Therefore, density 
fluctuations occur only at chain ends if no local moves are allowed. 

The main focus of this study is the dynamics of the "annealed SSA" (aSSA) where all chains are allowed to move: A 
chain and one of its ends are selected at random and a global move in a randomly chosen bond direction is attempted. 
We compare this with the "quenched SSA" (qSSA) of systems where all chains but one are quenched, thus, disallowing 
temporal density fluctuations (at chain ends or elsewhere). We sample the dynamics of the remaining free snake and 
average over all chains of one quenched equilibrium configuration. Unless stated otherwise, in the following we always 
refer to aSSA properties. Data from aSSA (qSSA) are gene rally indicated by open (full) symbols. 

Definition of a Monte Carlo step. Following Sokal jlj| we have chosen to set the unit of the Monte Carlo time 
step (MCS) such that on average each monomer is attempted once per cycle, irrespective of whether this attempt 
occurs through local or snake moves. If only the latter are allowed, every chain is on average visited once per MCS. 
Our definition of the time scale differs from the one used traditionally in the description of slithering snakes, which 
focuses on the computational time required to equilibrate samples. Indeed every snake move requires (essentially) the 
same CPU time irrespective of N, since all monomers are moved simultaneously. Hence, the computational relaxation 
time is proportional to r t /N where Tt stands for the terminal relaxation time in our units. 

III. PARAMETERS, CONFIGURATIONS, MEASUREMENTS 

This section presents various technical aspects of the simulation and introduces the important static and dynamical 
quantities for the subsequent analysis. When a quantity is defined, we often refer to the figures where its dependence 
on N, <)) or time is shown, without discussing the figure in detail, however. This discussion will be done in Sec. 

Parameters. In this paper we study athermal systems (no additional potential parameters) containing monodis- 
perse snakes (i.e., disallowing any local move) at volume fractions 4> ranging from dilute solutions to dense melts 
(4> ~ 0.5). The volume fraction is defined by <fr = 8NN p / L 3 , where N, N p and L are the chain length, the number of 
polymers in the system and the linear dimension of the simulation box, respectively. In the present work, we primarily 
investigate the scaling of dynamical properties with mass N (see, e.g., Figs. l2lfill9ll7ll8j) . Chain lengths ranging from 
N = 16 up to N = 1024 have been simulated. This is sufficient to put theoretical predictions to the test. We use a 
periodic cubic lattice of linear size L = 128. At (f) = 0.5 this corresponds to 131072 monomers per simulation box. 
Note that all length scales are given in units of the lattice spacing a (i.e., a = 1 in the following). 

In fact, we report here the first results from a larger and still on-going study, where in addition to N and 4> we have 
also varied the ratio to of local to snake moves over several orders of magnitude ranging continuously from purely local 
(lu = oo ) to pure snake dynamics (u) = 0). For the sake of comparison, the results of the local dynamics are included 
in some figures (Figs. I3I18|I . At the end of this paper, we briefly comment on the influence that a finite value for u 
has on the dynamics fFig.lTHJl. 

Measurement of static properties. Obviously, the static properties of an ergodic system do not depend on the 
dynamics of the algorithm. This has been carefully checked for the annealed SSA. (Numerically, it also appears to 
be true for the qSSA, even though ergodicity is less obvious in this case, since the environment of reference chains 
is frozen. Note however that all initial qSSA chains are taken from an equilibrated annealed ensemble.) This allows 
us to use the precise characterization obtained in previous studies, e.g., for the mean bond length I, the mean end- 
to-end distance R c or the mean radius of gyration R g . The first quantity gives the mean curvilinear distance along 
the chain associated with every successful snake move, the second the typical displacement R e /N of the center of 
mass per snake move. We recall that a polymer chain has a fractal dimension \ jv with v being the Flory exponent: 
R c oc R g ~ N v . The Flory exponent for dilute chains in a good solvent is v — vq sa 3/5. Above the overlap density 
<f>* 1=3 N/Rg 3 ~ 7V 1 ~ 3l/ °, the chains become Gaussian chains of blobs with v = 1/2. An important length scale for 
the subsequent discussion is the mean distance between chain ends, d G nd = (4-/V/</>) 1//3 , which we will often refer to as 
"mean end distance" in the following. Frequently, we use the "effective bond length" b{(j>) = lim,/v_voc R c /N u=1 / 2 
which is weakly density dependent. This definition is chain length independent only for systems above 4>* . In the 



dilute limit, we have R e /N u ° , which yields bo « 3 for N ^> 1. Some of these properties are summarized in Tables [I] and 
EH for (j) = 0.5 and = 0.125, respectively. Note that R g > d cnd = 27V 1 / 3 at <j> = 0.5 for N > 16, while for <f> = 0.125 
chains with N < 64 do not overlap, R g < d C nd- 

Various definitions of dynamic properties measured. As dynamical properties we have sampled various spatial and 
curvilinear mean-square displacements (MSD) gt(t), as well as associated times, diffusion coefficients and histograms. 
Following previous work |jj we define the spatial MSD g\{f) for the middle monomer (with monomer index n = N/2), 
92 (t) for the motion of the middle monomer in the frame of the center of mass, g% (t) for the center of mass motion, 
(?4 (t) for the end monomer motion (averaged over n — 1 or n — N) and 95 (t) for the motion of the end monomers 
in the center of mass frame. Examples are given in Figs. The MSD related to the curvilinear motion along the 
backbone is called ge(t) — (s(t) 2 ) where s is the curvilinear distance traveled by the middle monomer n = N/2 in 
time t. It is presented in Figs. 1121 1131 and H4l As we will see below, the scaling of ge(t) and the associated histograms 
(see Figs. I15I16I17|) are crucial for the understanding of the slithering snake dynamics. 

The curvilinear and spatial diffusion coefficients, presented in Tables HI and ITT1 and in Figs . Ifil 171 and IT81 are defined 
by D c = limt—xxj ge(t)/2t and D s = lim^oo g3(t)/6t. From the diffusion coefficients we introduce the time scales 
r s = R g 2 /D s and r c = (Nl) /(tt 2 D c ). The prefactor in the last definition has been introduced for convenience 
(see the discussion of Fig. El in Sec. IV Ej) . Additional times can be introduced to characterize the different MSD's. 
Following again Ref. the time scales n, T2, T3 and T4 are defined by g\(t = t\) = i? g 2 , g 2 (t = t 2 ) = 2i? g 2 /3, 
g 2 (t = T3) = g^(t = T3) and g§(t = t^) = R g 2 . In addition to this we shall use gi(t = T cn d) = d C nd 2 and g\ (t = r^) = d 2 
where d is an arbitrary distance. To characterize the curvilinear motion of the chains we also define the time t n / 2 f° r 
the diffusion of the middle monomer over half of the chain's contour by ge(t = t N / 2 ) = (Nl/2) 2 . From the analysis it 
will turn out that two time scales, ri and r on( j, are particularly important. They have been included in Tables [I] and 

m 

Finite size effects. Finally, we note that the results obtained for N — 1024 and (j> — 0.5 have to be taken with some 
care. The configuration contains only 128 chains and as L « R c ss 96, we cannot rule out dynamical finite size effects 
— especially for the quenched SSA. Some very large motions have been observed due to chains taking advantage of 
the image of their own correlation hole in the neighboring periodic boxes. Unfortunately, a systematic study of the 
influence of the box size L is beyond our computational capabilities at present. 



IV. ANALYTICAL EXPECTATIONS 



In the following paragraphes we summarize the relevant theoretical expectations for the analysis of the simulation 
data. 

Mean-square displacements: general properties. If the chains do not become permanently trapped, one expects free 
curvilinear and spatial diffusion, i.e., 

g 6 (t) ~ 2D c t (1) 
gi (t) =g 3 (t) = 9i {i) ~ 6D s t, (2) 

for times larger than the terminal relaxation time r t and for distances larger than the radius of gyration R g . Further- 
more, 32 and (75, both defined with respect to the center of mass of the chains, must yield a purely static quantity, 
the radius of gyration, 

sa(t) » Sb(*)/4 ~ R s 2 , (3) 

for very long times in any ergodic dynamical scheme @. 

At short times the spatial displacement is directly expressible in terms of the curvilinear displacement s along a 
fractal object characterized by the Flory exponent v. For the MSD of the center of mass this suggests 

1^1 )- { -^ 9 -f-m 2 ^) ( 4) 

where the penultimate step follows from a (seemlingly innocent) mean-field assumption. The last step uses the fact 
that in the melt v = 1/2. For the monomer displacements along the backbone a similar reasoning gives 

9i(t) - 92(t) « .94(0 - g 5 (t) « (R(s(t)) 2 ) « F b 2 (g 6 (t)/l 2 Y (5) 

where (i?(s) 2 ) is the size of a piece of chain of curvilinear length s. We have again made a mean-field approximation 
and assumed the fractal behavior of asymptotically long chains. 




Uncorrelated reptation hypothesis. Until now everything is fairly general and is valid for both the URH and the 
ARH. For uncorrelated snakes the motion along the backbone must be linear with respect to time. This fixes the short 
time predictions of the URH: 



Equation © refers, strictly speaking, to the approximately free diffusion along the chain backbone. As our simulation 
does not deal with phantom chains, there must be a small correction even for very dilute solutions. (An isolated chain 
is not a phantom chain. It still interacts with itself.) At long times the MSD's must obey Eqs. ^EJ). 

Activated reptation hypothesis for a quenched molecular field. The description of the MSD's within the ARH is 
best outlined by considering chains in a molecular field that is quenched at t = (qSSA). Hence, we start with 
an equilibrium ensemble of reference chains. This implies that the chains initially sit in their correlation hole with 
relatively low conformational free energies. To creep out of the local, well adapted environment a free energy penalty 
must be paid on average. It is then reasonable to assume (i) that fluctuations of the initial free energy may be 

neglected (i.e., we consider a typical initial chain), (ii) that the curvilinear displacement s may be taken as the reaction 
coordinate for the associated Kramers escape problem, and (iii) that the same penalty F(s) may be associated with the 
possibly different conformations characterized by the same s — a simple, but less obvious mean-field approximation. 

To estimate the dependence of the energy penalty on s it was argued in Refs. that F(s) = f s v for s < s + (N) 

where / is a constant. The upper limit s+{N) depends on the relaxation (escape) mechanism. Since we assumed the 
molecular field to be quenched, it seems natural to take s+ « TV ^> 1 : At distances larger than the radius of gyration, 
translational entropy must ultimately win over all other free energy contributions. To close the escape problem we 
still need to specify the physics for s > s+. Because of the gain of the translational entropy F(s) should decrease 
logarithmically at large distances. A crude, but physically sufficient approximation is to relax the penalty completely: 
F(s) = for s > s+ 0. 2 

The exponent ip appears to be a matter of debate. In Refs. 0,113 a linear relationship was suggested ((p = 1). 
This implies that the chemical potential for a curvilinear move is constant and thus independent of the tail length s. 
On the other hand, it was suggested in Ref. [23, that the free energy penalty should result from fluctuations of 
the molecular field rather than from a fixed chemical potential. The elaboration of this argument yields ip = 1/2. 

Given the free energy barrier one can calculate the mean-square displacements via Eqs. (1 lloll if one knows the prob- 
ability distribution p(s, t). In principle, p(s, t) can be obtained from the solution of the one-dimensional Smoluchowski 
equation. Analytically, this calculation has not been done yet, 3 but the qualitative behavior of the MSD's may be 
inferred from the following arguments. At short times [t < t_), we have s <C s+, and thus F(s)/k'oT <C 1. In this 
limit, one obviously recovers the free curvilinear diffusion of Eqs. © and l|7jl. If time increases, there should be an 
interval r_ <C t <C t + such that fcsT -C F(s) <C F(s+). Here, we expect that the chain is temporally localized in 
its correlation hole, which should lead to subdiffusive behavior due to intermittence of large scale displacements. The 
subdiffusive regime should extend up to a "hopping time" r + . r + is the time taken by the chain to diffuse over the 
curvilinear distance s+ and thus to overcome the activation barrier. This implies r + ~ s + 2 exp[F(s + )/k-QT]. Only for 
t^$> t + the chain can diffuse freely. 

Activated reptation hypothesis for an annealed molecular field. Up to now, we considered the molecular field to be 
quenched. If all chains are allowed to move simultaneously, this annealing should reduce the activation barrier. For 
this case it has been suggested that s+(N) ~ d en d 2 ~ iV 2 / 3 fTflL l2lll22T|: The free energy constraints can be relaxed 
by the reptation of other chains. For short times corresponding to s <C s+ the difference between a quenched or an 
annealed molecular field should be negligible. In this limit, the correlations between neighboring chains should be 
essentially of static nature. On the other hand, for s w s + the correlations are expected to be due to the cooperative 
rearrangements of several chains. Semenov and Rubinstein have proposed |22l | an intricate exchange mechanism which 
allows the relaxation of the free energy penalty between different snakes at s ss s+. It is not obvious to us if one can 
still cast these processes in the simple form of a one-dimensional Kramers escape problem. In any case this should 
only affect the probability distributions and MSD's around s ps s + . 

Relaxation times and diffusion coefficients. The relative scaling of the terminal time Tt and the diffusion coefficients 
is imposed by Eq. which fixes the ratio of both diffusion coefficients: 



g 6 {t)^2D c t, 
9\{t) = 92{t) ~ 9i{t) = g 5 (t) 



g 3 (t) « 6D s t 



(6) 
(7) 



(Nl) 2 _ Ej_ 
D c ~ D s 



(8) 



2 It is not sufficient to set F(s) = F(s+) for s > s+ as this corresponds to a much too high return probability for the one-dimensional 
Kramers problem into which we have cast the original three-dimensional model. 

3 We determined p(s,t) in the simulations, see Fig. 1161 



This ensures the matching of the long and short time behavior. Hence, it(iV), D C (N), D S (N) are expected to contain 
the same information. This should hold irrespective of whether the URH or the ARH applies. 4 

As we attempt to move all N monomers of a snake collectively, the curvilinear diffusion coefficient should be 
independent of chain length for uncorrelated snakes. Hence, it follows from Eq. @ and from the URH that 

N 2 

D c /l 2 Pa mN° , D B « m(R e /N) 2 , r t « , (9) 

m 

m being the N- independent mobility which, due to the URH, should be given by the acceptance rate, m pa A ~ N°. 
Note that these equations differ by a factor of N from the usual ones characteristic of reptating chains within the 
URH 3]. This is due to the fact that there, a factor N is spent by the internal Rouse dynamics along the primitive 
path. It is precisely this missing factor which makes the SSA attractive in the first place to explore the dynamics of 
extremely long chains. 

Within the ARH the curvilinear diffusion is fixed by the hopping time r+ to reach a distance s + within the potential 
formed by the local correlation hole, i.e., 

(s /I) 2 

D c (N)/l 2 re ' « mexp [- F(s+{N))] = mcxp( - fs+iN)*) = m+ , (10) 

which in turn implies D S (N) and it (A) according to Eq. (JSJ. Equation I|1U|) is obviously subject to the specific 
dependences of s+(N) and F(s) on their respective parameters. We will try to estimate both t+(N) and F(s) in 
Figs. [Hand CHI 

V. THE DYNAMICS OF ANNEALED AND QUENCHED SNAKE FIELDS 

We begin our presentation of the slithering snake dynamics in Sec. IV Al with the characterization of local phe- 
nomenological dynamical scales. Dilute and moderately dense systems are shown in Sec. IV Bl to be well described 
by the predictions of the URH. Diffusion coefficients and characteristic time scales for a broad range of densities are 
presented in Sec. IV Cl We then focus on the dynamical properties of systems of dense solutions and melts. The various 
MSD's and their scaling are discussed in the next subsection. By analyzing various displacement histograms possible 
ergodicity problems are explicitly ruled out for the annealed SSA in Sec. IV El The histograms of the curvilinear 
displacements allow us to estimate the free energy penality F(s) mentioned above. 

A. Local phenomenological dynamical parameters 

Acceptance rate. The simplest and statistically most accurate dynamical property sampled in any Monte Carlo 
simulation is the acceptance rate A (= ratio of accepted to proposed moves). It is included in the tables and plotted 
versus N in Fig. |3 and versus <f> in Fig. [3J We note en passant that the acceptance rates of the aSSA and the qSSA 
are identical. The SSA acceptance rates are also completely independent of chain length. We recall that in contrast 
a weak, but systematic chain length dependence of A is found for local dynamics due to the higher mobility of the 
chain ends compared to inner monomers (see, e.g., Fig. 10 of Ref. Q and the bold line in Fig. EJ. 

As shown in Fig. El ^ is density independent below (f> pa 0.01, where it is also much larger for the SSA than for 
purely local dynamics. This indicates that for these densities the blob size £ becomes sufficiently large to allow dilute 
SSA dynamics on distances smaller than £. Above <f> pa 0.1 density effects become pronounced and the acceptance 
rate of the SSA decreases strongly. This decrease may be explained as follows: one monomer is inserted in the SSA, 
i.e., 8 free lattice sites are needed, whereas 4 empty sites suffice for local moves, in which a monomer is shifted from 
its original position by one lattice constant. Hence, the observed smaller acceptance rate for the SSA at large density 
(0 = 0.5) is expected. 

Local mobility. Obviously, it would be naive to believe that the acceptance rate of the Monte Carlo process 
necessarily sets the scale for the local mobility of the monomers. There may be many allowed moves not contributing 



4 In Sec. IV Cl it will be checked explicitly that this general relation still holds for the annealed SSA at different densities (Fig.^J. As we 
shall see, it fails, however, for the quenched SSA at high density (Tigs, IT1I13I . 



to the motion. (They would also persist in a glass where the mobility is zero.) In order to measure directly the 
monomeric mobility m in the SSA we exploit Eq. Q and compare the simulation data with 

gi (t) = b 2 (mty. (11) 

This equation provides a reasonable description of the MSD of the central monomer at short times, as Fig. ^illustrates 
for a dilute solution and Fig. for a dense melt. 

For dilute systems it is appropriate to fit gi(t) using v = vq and the density independent effective bond length 
&o ~ 3. The same analysis can also be done in semidilute solutions for times shorter than the time required to cross 
the blob size so that the excluded volume interaction is not yet screened. This approach was also applied in Ref. Q 
to extract m for the BFM in the case when only local moves are allowed. Our results for m for both the annealed and 
the quenched SSA are included in Fig. El 

However, Eq. (|llf) with b = bo and v = vq is not the only possible way to do the analysis. For overlapping systems 
another choice is b — b(<p) and v = 1/2. This measures the effective mobility on distances of several blob sizes. 
Equation (|llf> then suggests that gi(t) should scale as gi(t) ~ b{4>) 2 t 1 / 2 at short times, independent of chain length. 
Figure El confirms this expectation. If the mobility is now extracted from Fig. EI the differences between the new values 
and those obtained before (using b = bo and v — vo) are negligible. 

Here, the important point is that the local mobility is the same for all chain lengths at a given density. (Note that 
the mobilities obtained by Paul et al. 7} for the local dynamics are also independent of chain length.) Furthermore, 
Figs. El and El also show that m is independent of whether the system is annealed (aSSA) or quenched (qSSA). For 
both dynamical schemes, the mobility rapidly decreases with density if cf> > 0.1. This behavior is qualitatively the 
same as that of the acceptance rate. Similar results are obtained if Eq. Qlljl is fitted to the MSD of the chain ends 
34 (i) (not shown). For the SSA, the chain ends are more mobile than the central monomer by a factor of order one. 

In summary, the acceptance rate or the local mobility cannot explain the unexpected chain length dependence of 
the motion of strongly entangled chains to be discussed in Sees. IV CHV El 

Effective mobility for the aSSA. A glance at Fig. El shows that in fact two fits using Eq. (jl 1|> are possible for 
the annealed SSA once <fi ^ </>*■ The rescaling y = gi(t) /b(<j)) 2 t v=1 l 2 yields a data collapse at very short times, 
t -C t_. From the resulting plateau the local mobility m, shown in Fig. El was extracted. For larger times, however, y 
decreases within the window r_<t<C t + (N). The existence of such a time window is obviously at variance with the 
uncorrelated reptation hypothesis. For still larger times, r + -C t <C r t , the rescaled MSD becomes approximatively 
horizontal. This allows us to define an effective mobility m + (N) (included in the tables). It reflects the complicated 
correlations in the second regime r_ <C t <C r+, which are responsible for its pronounced dependence on N. The 
importance of m+ for the description of diffusion coefficients and relaxation times will be clarified in Sec. lVCl Finally 
for t > r t , Fig. El shows that y cx yt, implying that the monomers diffuse freely. 

B. Dilute and weakly overlapping systems 

Before analyzing the intricate physics at high chain overlap we briefly consider the simpler situation of dilute or 
weakly overlapping solutions. These systems serve as a reference point for the following discussion of concentrated 
solutions and melts. That the predictions of the URH apply in the limit of high dilution can be inferred from Figs. El 
and El The rescaling used for the vertical axes of the three figures is motivated by Eq. © and the fact that the local 
mobility m w A(d>) ~ N°. The figures confirm the validity of the URH by showing that D c ~ N°, ND S ~ N 2 " ' 1 
and r t « n ~ iV 2 . 

Figure El presents the spatial MSD's for a typical dilute or weakly overlapping configuration. We compare the 
measured MSD's with the expected long time behavior (Eqs. I|2I3|I ) and the short time predictions of the URH 
(Eqs. (|6I7|I using v = vq). Note that the MSD of the center of mass 33 (t) is linear over several orders of magnitude, 
implying that corrections due to excluded volume correlations of the snake with itself are very small. The same 
applies to the curvilinear motion gg (t) (not shown). Incidentally, the excellent agreement confirms that our program 
is running properly and that the URH is appropriate as long as chain interactions are weak. The presented data are 
for the annealed SSA, but absolutely identical results have been obtained for qSSA, as expected. 

The above statements remain valid for all configurations and both dynamical schemes provided that (j> < <f>*. 

C. Diffusion constants and time scales 

Diffusion coefficients for the aSSA. The curvilinear and spatial diffusion coefficients, D c and ND S , are displayed 
in Figs. El and El not only for dilute solutions, but also for densities ranging up to the melt density <j> — 0.5. We plotted 
D c /A and D S N/Ab 2 versus chain length, since we expected to find more or less horizontal curves (see Eq. |(jJJ). 



Surprisingly, Eq. @ breaks down with increasing chain overlap. For a given density the data points become pro- 
gressively curved in log-log coordinates. In fact, D c and ND S drop over more than two decades at <fi = 0.5, as N 
increases. It is important to stress that the chain length dependence of D c and ND S increases gradually with chain 
overlap and that this effect is by no means a pathology of a lattice model at high density. 

The data for cf> = 0.5 may be compared to a phenomenological power law with exponent j3 — 4/3. This power law is 
merely a guide to the eye and a rough attempt to characterize the slope in log-log coordinates for the available chain 
lengths. However, the stretched exponential D c sw ND S w cxp(— 0.87V 1 / 3 ) motivated by Eq. (| 1 Of) describes the data 
much better. Our scaling is, in fact, in excellent agreement with the ARH and the proposed energy barrier exponent 
ip = 1/2 [2]]] • The alternatively suggested value tp = 1 0,H2 is not compatible with the simulation data (at least for 
the chain lengths simulated up to now). 

If, instead of the acceptance rate, we use m + {N) to rescale the diffusion coefficient, i.e., ND s /m+b((j)) 2 , we find 
a very good data collapse (small symbols at the top of Fig. Q). A similar collapse is also possible for D c /m + (not 
shown). This shows that it is possible to absorb the complicated physics at short times into one parameter, m + (N). 
Thus, the dynamics for r + <C t -C r t is perfectly described in terms of uncorrelated reptation dynamics along the 
fractal chains. It remains to describe and possibly explain the dynamics in the second regime, r_ <C t -C r + , of Fig.0 

Diffusion coefficients for the qSSA. With increasing chain overlap the qSSA diffusion coefficients decrease even 
more rapidly with density and chain length than their aSSA counterparts. In fact, while we have obtained most of 
the curvilinear diffusion coefficients from our qSSA simulations, we have been unable to measure D s for (f> > 0.125. 
This surprising statement is obviously at variance with Eq. I[8"|l. It indicates that curvilinear and spatial motions 
decouple for these high densities due to the localization of the snakes within their correlation hole in the quenched 
case. Furthermore, at high density D c is found to be a power law without any detectable curvature (Fig.EJ). However, 
it is interesting to check whether Eq. (|1C)|> with s+(N) = N for the qSSA and p = 1/2 (for both schemes) yields a 
reasonable fit. This is indeed the case: D c /A ss 20 exp(— 0.15N 1 / 2 ) is hardly distinguishable from the indicated slope 
P 1.5. Longer chains are obviously necessary to confirm the exponential behavior expected by the ARH. 

Relaxation times for the aSSA. In Figs.|Hland|niwe display some of the characteristic times defined at the end of 
Sec. IIIII In the first figure we test Eq. JHJ by plotting various ratios tJt\. Apparently, the ratios are independent of 
chain length. We find r s /4 as r c /5 wti« T2/O.85 « T3/4.7 T4/O.3 t^/ 2 /&. Only the data for <fi = 0.5 have been 
included, but the same results prevail for all densities. This demonstrates that all indicated times are proportional 
to one characteristic time scale, the "terminal time" r t , containing the same information as the diffusion coefficients. 
That is, Eq. © still holds, although Eq. © fails. 

Accordingly, Fig. |^1 shows, e.g., that t\ is proportional to N 2 in the dilute limit, but increases like 7*1 « N 2 /D c ~ 
iV 2 exp(0.87V 1 / 3 ) in the melt. The fi gure also includes for d = 3 and r en d- The first time indicates the distance d 
up to which c/i (t) remains independent of N, and characterizes the extent of the first dynamical regime displayed in 
Figs. IB1 and 1121 In fact, it strongly decreases with monomer density, and we expect it to scale like the time needed 
to diffuse over the excluded volume blob £. From the observed behavior of D c one anticipates for r cn d the scaling 
n ^> T cn d ~ s cnd/^c ~ N 2 ^ 11 exp(0.8A^ 1 / 3 ), where s cn d ~ rfond 1 ^- This is again consistent with the data presented in 
the figure. 

Relaxations times for the qSSA. Data points for t% and t n / 2 at <j> = 0.5 have been included for the qSSA in Fig. [5] 
Unfortunately, the range of N, which could be simulated, is insufficient to check if indeed t\ ~ N 2 exp(constA^ 1 / 2 ). 
The time scale tjv/2, measuring the curvilinear motion over N/2 monomers, shows an effective power law which is 
consistent with the corresponding diffusion coefficient from Fig. El Tyv/2 ~ N 2 /D c ~ iV 2+15 . Note that T\ and T7V/2 
scale differently in contrast to what was observed for the annealed SSA (Fig. EJ. 

D. Mean-square displacements of dense snakes 

Long time limit. Focusing on systems in the melt limit (<f> = 0.5) we now discuss several MSD's and their scaling 
in detail (Figs. As a glance at the presented MSD's shows, the long time predictions (Eq. (^OJ) are in very 

good agreement with the data for the aSSA, indicating that the chains do not get trapped and that the simulation 
runs have been long enough. This can be seen in Fig.^Jfor the spatial MSD's for N = 512, in Fig.^Jfor both spatial 
and curvilinear MSD's for N = 128 and in Fig. U3 for the curvilinear MSD g 6 (t). 

The situation is more complex for the data of the qSSA (included in Figs. llllll2J) . Figure ITTI shows that the spatial 
MSD's increase (approximately) logarithmically for the longest times we have been able to simulate. This must be a 
transient: For very long times these MSD's must either become constant, i.e., the chains become localized on a length 
scale (presumably) given by the radius of gyration, or ultimately, they diffuse freely in space. So far, we have been 
able to show conclusively only the correctness of the second proposition for small chains (TV < 32) . 

At first sight surprisingly, it is much less difficult to reach the free curvilinear diffusion limit, as may be seen from 
Fig- El where ge(t) and Ngs(t) are directly compared for N — 128, and from Fig. ^1 where the g§(t) for different 



masses are presented. This indicates that for the qSSA the curvilinear and spatial MSD's decouple for t 3> T coup and 
s S> s C oup (see Fig. lll(l : although unable to move much further in space, the snakes manage to move relatively rapidly 
within the cage of their correlation hole. This is why we prefer to call this surprising behavior of dense qSSA snakes 
(permanent or transient) localization in contrast to a trivial local blocking of the snake moves (on a scale given by the 
monomer density), which must obviously occur at very high densities. 

Short time limit for the aSSA. While the MSD's are well described by the power laws of the URH in the dilute 
limit (Fig. IU, this becomes strikingly different at larger densities in the time window t_ < ( < t + . This is revealed 
by the direct comparison of dilute and dense systems with N = 512, as given in Figs. l4l and ITU1 In this paper we 
attempt to characterize and to explain the features visible in the latter plot. 

We have tentatively described the subdiffusive behavior of the aSSA data by means of an effective power law: 



This power law provides a remarkable fit for r_ <C t <C t + (compare Figs. l5Tl0ll2(l . even though some curvature is 
evident in log-log coordinates. The alternative fit 



is less satisfactory for our aSSA data (fitting a smaller time interval), as may be seen from Fig. ^2 (The exponents 
in the two previous equations have naturally been chosen consistently with the general relationship expressed by 
Eqs. |g|3 using v = 1/2.) 

Short time limit for the qSSA. Equation (IB) , however, reasonably fits the qSSA data. This may be seen from 
Figs. IB1 and 1111 While for the spatial MSD's this fit works for all times t ^S> t_, ge(t) becomes linear for t 3> r +q , 
t 3> r +a with r +q 3> r +a (compare aSSA and qSSA in Fig. ^]for N = 64, 128 and 1024), where we altered the index 
to distinguish the different algorithms. Note that annealed SSA and quenched SSA data are identical for small times 
t <C Tfl uc (Af) (Tigs. ITlll2(l . This is due to the low probability that a snake in the aSSA may interact with the density 
fluctuations created by neighboring snakes at short times. Hence, the molecular fields probed behave in the same way. 
As Tfl uc (N) increases strongly with N (see below) , this suggests that the logarithmic behavior might be seen for aSSA 
snakes as well, if we could simulate longer chains. 

Scaling of g\(t), g 3 {t) versus ge(t). In Fig. 1181 we verify now directly whether the time dependence of all spatial 
MSD's is indeed fully encapsulated in terms of the curvilinear MSD ge(t), as suggested by both the URH and the ARH. 
In the main panel of the figure we check this for the aSSA data at <f> = 0.5 by plotting gi(t) and ga(t) versus ge(t). 
The MSD of the center of mass is linear with regard to g§, as expected. This is in accord with the previous finding 
D c w ND S for Gaussian chains. We find a similar data collapse for gi(ge) for i < T t . For sufficiently long chains, we 
clearly recover the anticipated power law of Eq. (7J for the diffusion on an object of inverse fractal dimension v = 1/2 
(dashed line in Fig. 1 1 31) . (The departure from the slope 1/2 for small g&(t) is due to the non-gaussian statistics for 
small distances along the chain.) We conclude that the curvilinear and spatial diffusion are translated into each other, 
as expected, and that the departure from the URH can be traced back to the anomalous curvilinear diffusion depicted 
in Fig. [H 

In the inset of Fig. ^| we present Ng 3 (g§) for the qSSA for two chain lengths and two densities. While for the 
lower density g 3 {t) remains a (roughly) linear function of ge(t), this is only true for short times at cj) = 0.5, where 
the spatial displacements increase much more strongly. Note that both densities appear to be independent of chain 
length. This indicates a very weak mass dependence for s COU p- 

Scaling with chain length. We now discuss in more detail the chain length dependence of the subdiffusive behavior 
at fixed volume fraction <fi = 0.5. Since we have shown (Tig. 113(1 that the spatial MSD's of the annealed SSA scale in 
terms ol ge(t), we can concentrate on the N dependence of the latter quantity (Tig. ll2[l . The vertical axis is normalized 
to yield the curvilinear diffusion coefficient D c as the plateau value for long times. At short times t <C r_ w 50, all 
curves merge. Hence, r_ is chain length independent, but depends weakly on the density, as may be shown by a similar 
plot containing data for different </>. For times beyond r_ the dynamics slows down more and more until it reaches the 
plateau of free diffusion at r +a (A) for aSSA and T +q (A) for qSSA, respectively. Apparently, the correlations become 
more and more pronounced with increasing N. 

We have directly computed the four times r +a , r +q , r nuc and r coup and the related curvilinear displacements. Our 
measurements are not very precise (particularly, for the qSSA data better statistics is needed), and we only present 
some general trends here. A comparison of all chain lengths shows that 



Ng 3 (t) w 9l (t) 2 w 6 (t) ~ t°- with awl/2. 



(12) 



Ng 3 (t) ^ 9l (tf « 56 (i) ~ log(t) 



(13) 



Tfluc (AO 
Tcoup (A) 



T+ a (A) « T cnd (A) 
T +q (A) ~ N 2 . 



(14) 

(15) 



In other words, both pairs of time scales correspond to different definitions of the same physical processes. In the first 
case, these are the density fluctuations in the annealed scheme, causing aSSA and qSSA to differ and the curvilinear 



motion to become eventually diffusive. The second equation corresponds to the free curvilinear diffusion of the snake 
trapped in its quenched correlation hole. The power law indicated in Eq. H15[l is an estimate motivated by the data 
given in Fig. El (diamonds). 

The last claim of Eq. 114|) , which is central to the understanding of the observed correlations and is in line with 
the ARH [T^. |2^ | , is further tested in Fig. The attempt r +a w t\ fails to scale the crossover regime between 
the subdiffusive time window and the free curvilinear diffusion at late times. (Note that there is always the time 
scale t_ ~ N Q for the free diffusion at very short times and no scaling covering all regimes is possible.) Instead, the 
alternative r +a w r cn d is successful. The observed scaling confirms Eq. (|14H and demonstrates that the crossover from 
correlated to mean- field behavior is determined by d en d- 

Finally, we note for the second, less important, relaxation process for the qSSA that s +q increases very weakly 
with mass, s^ q ~ A - 25 . (However, this does not contradict the scaling documented in the inset of Fig. Note 
also that this is again a rough estimate, but we definitely have to rule out the two "natural" guesses s+ q ~ N or 
s +q ~ N°.) Interestingly, s^ q ~ TV 0,25 and r +q (iV) ~ TV 2 are compatible with the (very precisely measured) qSSA 
diffusion coefficient from Fig.® D c s^ q /r +q ~ 1/7V 15 , and, hence, with the time scale ~ A 3 ' 5 shown in Fig.0 
While this is satisfactory, we presently cannot explain the scaling behavior of r +q and s+q. 

E. Histograms and Non-Gaussianity 

So far, we have discussed results obtained from the sampled mean-square displacements and related diffusion 
constants and time scales. We now turn to the additional information contained in the histograms and higher moments 
of the curvilinear displacements. Merely results from the annealed field dynamics are discussed here, since statistics 
is too poor for the qSSA. 

An important, not only technical, question concerns the possibility that some chains or monomers become locally 
completely blocked. We have paid much attention to this issue and have carefully checked that, in fact, all chains and 
monomers ultimately move arbitrarily far in space for the densities studied (<p < 0.5). Clearly, for some higher density 
this must eventually break down. 

One may readily construct by hand chain configurations of the BFM, which lead to blocked chains for any dynamical 
scheme under consideration, even if local moves are allowed. These extremely rare points of the phase space also exist 
for our slithering snake scheme, although they are apparently excluded from our initial configurations (which is 
the physically important point to make here). These states cannot be accessed (by construction) and in this strict 
mathematical sense our dynamics is not ergodic. However, this is completely irrelevant for the thermodynamic and 
dynamic behavior. 5 

Histogram of first passage time. Two checks are presented in the histograms given in Figs. 1151 and 1161 The first 
figure presents the histograms of first passage times of the central monomer over a curvilinear distance s = IN/2 
for various chain lengths and densities. At long enough times, typically 10rf p , where Tf p = (if p ) is the mean first 
passage time, all chains have crept out of their initial tube. If properly rescaled, the data collapse onto a master curve. 
The distribution of the first passage times can be calculated if free curvilinear diffusion is assumed. It is a linear 
superposition of modes p with an exponential time dependence exp(— p 2 t[ p /T c ). For long times, this sum is dominated 
by the lowest mode and, accordingly, we find an exponential decrease exp(— if p /r c ) for all systems. As expected, the 
mean first passage time scales like the terminal time, Tf p w 6ti « 5r c /6, and is of the same order as t n / 2 which was 
defined as the time where gs(t) corresponds to the same curvilinear distance. 

Histogram of curvilinear displacements s. Figure ITfil presents the histogram of s(At) for N — 1024 and <f> — 0.5 for 
various time intervals, as indicated. We focus on the interesting time window r_ «KT + «T t . First of all we stress 
that the histograms have (very qualitatively) the same shape and fall off monotonously without any singularity at 
the origin or elsewhere, i.e., no snake gets blocked at small s. If no correlations were present, the rescaled histograms 
should all collapse onto a single Gaussian. More importantly, however, Fig. ^5] shows that the different dynamical 
regimes correspond to (slightly) different distributions. While at times t <§C t_ and t t + (N — 1024) w 0.2 ■ 10 8 
the distributions are in fact Gaussian 6 , they are essentially exponential in the subdiffusive regime where density 
fluctuations are important. The distribution is thus much broader than originally expected, but all moments of the 
distribution exist. Although this is not an equilibrium distribution measuring directly the free energy as a function of 



In this context it is also of relevance that one obtains the same dynamics if a small amount of local hopping moves, u) <C 1, is included 
in the algorithm which allows additional transitions between points of the phase space. See Fig. 1181 

For clarity, both regimes are not included in the figure. Note that Gaussian behavior is only found in the first regime at much smaller 
densities than the one discussed in Fig. 1161 A sufficiently large blob is required to obtain a Gaussian distribution. 



s, the Smoluchowski equation suggests that — log(p(s)) should be proportional F(s). We have therefore compared the 
two proposed exponents for the activation penality ip — 1 (bold line) and ip = 1/2 (dots) with the data. Unfortunately 
at large s, the statistics becomes too poor to distinguish between both exponents. More work is required to relate the 
form and scaling of the observed histograms to the respective MSD. 

Non-Gaussian parameter. The non-gaussian parameter a 2 ±$, - 1 gives a compact means to describe the 
Gaussianity of these distributions as a function of time. This is presented in Fig. El We recall that a 2 — for a 
Gaussian and a 2 = 1 for an exponential distribution. At short times a 2 (t) decreases rapidly and is chain length 
independent. The non-gaussianity at these times is due to local monomer interactions. In agreement with the his- 
tograms shown in Fig. 1161 a 2 (t) increases strongly for larger times. This shows again that correlations become more 
pronounced. The effect increases strongly with chain length in agreement with the subdiffusive behavior in Fig. 1121 
Note that a 2 becomes even much larger than 1 for our largest chain N = 1024, reflecting the existence of a broad tail 
in the distribution shown in Fig. 1161 We caution, however, that finite size effects might contribute to this effect, as 
the end-to-end distance of this configuration is comparable to the box-size. It is possible that some snakes manage to 
find their way rapidly through the periodic box by taking advantage of their own groove spanning the box. For times 
larger than t+(N) the non-gaussianity vanishes again for all N, as it should. 

VI. SUMMARY AND OUTLOOK: SNAKE VERSUS LOCAL MOVES 

Summary. In this paper we reported first numerical results on the dynamics of monodispcrse polymer chains 
moving according to the slithering snake algorithm SSA (no local relaxation pathways). The simulations were done 
with the three-dimensional bond- fluctuation model (BFM). We investigated two different dynamical schemes by 
focusing mainly on the polymer melts (high monomer densities): Either we allow all snakes to move (aSSA), which 
permits fluctuations of the molecular field at the ends of neighboring chains, or all snakes but one are quenched (qSSA). 
This work concentrates on the annealed algorithm, whereas the second mainly serves as a benchmark to investigate 
the role of density fluctuations due to neighboring chains. We studied the scaling of various spatial and curvilinear 
mean-square displacements (MSD), diffusion coefficients and relaxation times with chain length. Furthermore, we also 
presented the distributions of the curvilinear motion. 

Our results are broadly in agreement with the "activated reptation" hypothesis (ARH) |TgLl2^. For the aSSA, we 
find that at high chain overlap the terminal relaxation time Tt increases much more strongly with chain length N 
than possible for uncorrelated snakes which are obtained for dilute and weakly overlapping chains. In fact, a stretched 
exponential r t ~ N 2 exp(0.8iV 1 / 3 ) fits our data well. Similar and directly related (Eq. JHJl) behavior can be found 
for the spatial and curvilinear diffusion coefficients, ND S ss D c , and other time scales, such as r cn( j, which measures 
the time needed for the diffusion over the mean distance between chain ends d on d- The stretching exponent 1/3 is in 
agreement with the free energy penalty F(s) ~ s 1 / 2 (s being the curvilinear displacement) estimated in Ref. [21|, but 
is incompatible with the more recent suggestion F(s) ~ s of Ref. [2^ ( a t least for the chain lengths accessible in the 
present simulation). 

A more detailed study of different moments of the displacements shows that the snake dynamics is only uncorrelated 
(exhibiting a free curvilinear diffusion with Gaussian distribution) for times smaller than T-(<f>) ~ N°, which is 
determined by the blob size £, and for times larger than t + w r en( j (shown by the scaling collaps of Fig. I14|l . which 
is due to the density fluctuations present for distances larger than d en d- They are responsible for an additional time 
scale (other than the terminal time) in the annealed algorithm, allowing the chains to relax the free energy penality 
caused by their motion in a locally rigid and frozen network of obstacles retaining them in their correlation hole. This 
is supported by the finding that the MSD's of both the aSSA and the qSSA, which are strictly identical at small 
times, become different around r + . 

In the intermediate time interval the dynamics becomes increasingly subdiffusive and non-gaussian fFigs. ITf)ll7|l . 
as the chain length is increased for both aSSA and qSSA dynamics. The precise description of the MSD's in this 
regime is not evident for the masses we have been able to simulate and our study is not fully conclusive. While power 
laws yield nominally better fits for aSSA, logarithmically slow relaxation processes might also be appropriate. The 
observed behavior for the quenched dynamics fFigs. ITlll2(l . which could provide an "envelope" that the annealed 
dynamics might follow more closely with increasing chain length, is indeed better described in this way. Logarithmic 
behavior is, however, not easily reconciled with (mean-field) ARH approach which predicts an algebraic relationship. 
More simulations are warranted to clarify this issue. 

The curvilinear displacement histograms decrease monotonously without any singularity at small s. In the subdif- 
fusive intermediate time they decrease essentially exponentially (Fig. Ilfijl. This suggests that the free energy penalty 
F(s) is indeed of power law form. A direct measurement of the acceptance rate to perform a step s — > s + 1 is currently 
underway using the quenched scheme which should allow a direct computation of the free energy cost for an escape 
attempt of order s. 



In summary, we believe that our study demonstrates that the activated reptation approach is relevant for the 
description of concentrated solutions and melts of purely slithering snakes without local motion (ui = Q). 

The mapping onto a locally realistic model. It is important to recall that the slithering snake dynamics discussed 
in this paper is highly artificial, as no local rearrangements are allowed and density fluctuations can only occur at 
chain ends. The idea is that the snake corresponds to the primitive path of the reptation model and each monomer to 
a coarse-grained tube segment of an originally local model. The mapping of both levels of description is non-trivial. 
It is not evident whether it can be performed at all without an additional tuning parameter, such as the frequency of 
local hopping moves. 

The mapping might involve the density. It is possible that some of the high molecular densities, which we studied 
here, are in fact too large to correspond to realistic volume densities of an underlying microscopic model. A character- 
istic density might therefore exist below which even asymptotically long snakes remain uncorrelated. At present, our 
computational results do not seem to favor this option, as suggested, e.g., by Fig. [7J However, it remains a possibility 
which is currently pursued using lower densities and much larger chains such that 1 > ^ > cb*(N). In any case, a 
more detailed analysis of the density dependence is warranted to establish the scaling of r_ (0) and t+ (N, </>) . In other 
words, it has to be confirmed unambiguously that the first time scale is related to the blob size £((/>) and the second 
one to the mean distance between chain ends d en d(N, 4>). 

On the importance of local motion. Another issue is related to the fact that our slithering snakes do not have any 
mechanism to relax constraints via local pathways. It might well be that some local dynamical degrees of freedom are 
necessary to tune a mesoscopic computation scheme such that it becomes consistent with real reptational dynamics. 
A finite fraction of local moves should weaken the strength of the obstacles imposed by neighboring chains. This 
is indeed borne out by preliminary studies presented in Fig. 1181 where the rescaled spatial diffusion coefficient is 
plotted for various ratios of local to snake moves ui. For short chains (N < 64) ND S decreases monotonously with 
increasing u, since local moves are less efficient in exploring the phase space and the confinement is neglible. As 
the mass increases, we find a striking non-monotonous behavior. The dynamics first becomes more rapid, as local 
dense fluctuations become more relevant. This effect apparently saturates at a frequency of order 10, which causes 
the diffusion coefficient to decrease again strongly with increasing u> (at fixed N). A finite ui generates a lateral tube 
of size d\(ui, </>) > a, which should become identical to the tube diameter d e (4>) of the reptation model in the limit of 
very large u>. Recalling that d c (cf>) s» £(<fi) we speculate that a successful mapping of a SSA scheme on a local scheme 
may require the tuning of the ratio di(u), (f))/£(4>) — const « 1, which implies a different u> for each density. We are 
currently pursuing further simulations to clarify these questions. 

A final word. It might well turn out that even in this larger modeling space no lv(4>) can be found, where D c ss ND S 
remains chain length independent for large N. We would regard this as a strong evidence for the correctness of the 
ARH for a model with fully realistic local dynamics. We believe, however, that astronomically long chains are required 
to obtain measurable effects. 
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I 


b(N) 


A 


m+ 


ND S 


D c 


n 


Tend 


16 


2.602 


2.9 


0.08433 


0.0164 


0.057 


0.136 


299 


376 


32 


2.603 


3.0 


0.08478 


0.0123 


0.040 


0.085 


1588 


1257 


64 


2.603 


3.1 


0.08502 


0.0071 


0.022 


0.045 


14527 


5434 


128 


2.603 


3.1 


0.08515 


0.0035 


0.010 


0.021 


1.3 ■ 10 5 


27826 


256 


2.603 


3.2 


0.08506 


0.0012 


0.0037 


0.0075 


1.6 ■ 10 6 


192900 


512 


2.602 


3.2 


0.08637 


0.0003 


0.0010 


0.0021 


1.8 • 10 7 


0.4 • 10 7 


1024 


2.603 


3.0 


0.08606 


0.0001 


0.0002 


0.0005 


2.5 • 10 s 


0.2 • 10 s 



TABLE I: Various static and dynamical properties for slithering snakes with aSSA dynamics at volume fraction 4> = 0.5 versus 
chain length N: The mean bond length I, the effective bond length b(N) = Rc/N 1 ^ 2 , the acceptance rate A (Fig. the 
effective mobility m+(N) from gi(t) (Fig. 0, the spatial self-diffusion constant D s = limt^oo gs(t)/6t (Fig. 0, the curvilinear 
diffusion constant D c — limt_ 00 ge(t)/2t (Fig. and, finally, the time scales n and r cn d obtained from <7i(ti) = i? g 2 and 
ffi(T"cnd) = 4nd 2 , respectively (Fig. [5J. The first of these times indicates the terminal time (ri ~ r t ), the latter separates the 
strongly correlated anomalous chain diffusion from the free curvilinear diffusion at longer times (r cn d ~ r+). The mean bond 
length / and the acceptance rate A are essentially chain length independent while the effective mobility m+ clearly is not. 
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116 
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0.62493 
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0.86 
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94 


188 
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2.691 


3.7 


0.62759 


0.13 


0.79 


1.2 


533 


743 


128 


2.692 


3.8 


0.62888 


0.11 


0.66 


1.0 


2246 


1680 


256 


2.692 


3.9 


0.62940 


0.09 


0.57 


0.8 


15142 


6425 


512 


2.692 


3.9 


0.63008 


0.085 


0.36 


0.5 


94401 


20059 


1024 


2.692 


4.0 


0.63006 


0.055 


0.28 


0.3 


780000 


26329 



TABLE II: Same properties as in Table |T] but for a lower volume fraction cj> — 0.125. Note that the effective bond length 
b(N) = Rc/N 1 / 2 increases slightly with chain length: The chain overlap is weaker than at <j> = 0.5 and larger chains are needed 
to screen the excluded volume correlations. For the relatively dilute systems < 64 the time scale r cn d becomes larger than 
the terminal relaxation time n. 




FIG. 1: Sketch of the bond-fluctuation model (BFM) with local (top chain) and slithering snake moves (bottom chain). We 
show a two dimensional projection of the three dimensional Monte Carlo algorithm. Monomers (squares) are connected by 108 
possible bonds (bold lines) allowing for 87 angles between bonds (in three dimensions). In the original BFM, monomers make 
only local jumps in six spatial directions, in (automatic) agreement with the imposed topological constraints. In this paper we 
study the dynamics of the slithering snake algorithm (SSA). As depicted for the bottom chain, monomers move collectively 
along the chain in the SSA. Effectively, this amounts to removing a monomer (the striped one on the left), connecting it to 
the other end of the chain and leaving the middle monomers unchanged. Therefore, density fluctuations and constraint release 
only occur at the chain ends. 
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FIG. 2: Acceptance rate A versus N for various densities <j>, as indicated in the figure. The large open symbols correspond to 
the annealed SSA dynamics, while the small, filled symbols indicate the quenched SSA data. Both results are identical. For the 
SSA the acceptance rate is independent of chain length. Contrary to that, the bold line indicates A at <j> = 0.5 for the BFM 
with purely local moves. Here, a slight dependence on N is found. 
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FIG. 3: Acceptance rate A (spheres) and local mobility m (squares) versus <f> for pure slithering snake moves (u> = 0): annealed 
SSA (open symbols) and quenched SSA (small, filled symbols). Both A and m are not affected by the freezing of the molecular 
field. Also included are the acceptance rates from the BFM with local moves only (u> — oo), extrapolated to infinite chain 
length: limj\r-»oc A(N, (j>) (bold line). For the SSA the acceptance rate and the local mobility are independent of N. 
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FIG. 4: Different spatial mean-square displacements (MSD) g%(t) for a dilute system of chain length N = 512 at volume fraction 
cj> — 0.0156. The MSD's are defined on page0] Only the data for the annealed SSA are given, as both dynamical schemes (aSSA 
and qSSA) yield the same results in the dilute limit. At short times, gi(t) = gift) ~ gi(t) = gs(t) oc ^"= - 588 (bold solid lines), 
as expected from Eq. £J. In good agreement with Eq. we also find gz{t) oc t/N over (almost) the entire time range. The 
straight (dashed) line for gz(t) is a fit to the long time asymptotic limit. For large times gi(t) and g4,(t) must merge with gz (t), 
while g2(t) — > (Rg) and gs(t) — > i(Rg) converge towards time independent static properties. 
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FIG. 5: Estimation of local and effective mobilities, m and m+. The rescaled MSD y = 31 (t)/(fo(0) 2 t 1,/2 ) for the displacements 
of the central monomer is shown at the melt density cj> — 0.5. The mass N varies, as indicated. Open symbols refer to the 
annealed SSA. We also included the qSSA results for N = 64 and N = 128 (filled symbols). For comparison, we show a power 
law with exponent —1/4 (dashed line) and gi(t) « 3.51og(t) — 3.8 (bold line) which fits well the subdiffusive behavior of the 
qSSA. From the two horizontal envelopes of the curves we estimate the short time mobility rn((j>) ~ N° and the effective 
mobility of the snake dynamics m+(N, <f>). (The latter quantity is not visible for qSSA, as a much longer simulation run would 
have been required.) As can be seen from the figure, m+ strongly depends on chain length. 




FIG. 6: Curvilinear diffusion coefficient D c /A versus N for various volume fractions <f> (A = acceptance rate). The open 
symbols show the results for the annealed SSA, the filled symbols for the quenched SSA. Only in the dilute limit D c /A ~ iV , 
as expected. Although the effect strongly increases with <j>, we find deviations from this expectation at all densities for large 
chain overlap <j>/(f>*(N) 3> 1. The effective power law with exponent f3 = 4/3 roughly characterizes the chain length dependence 
at tf> = 0.5 for the aSSA dynamics. A better fit, which properly takes into account the curvature of the data, is given by the 
stretched exponential exp(— 0.8A" 1//3 ) (bold line). The qSSA data at <f> = 0.5 are well fitted by an effective power law (/3 ~ 1.5) 
without any noticeable curvature (dashed-dotted line). 




FIG. 7: Spatial diffusion coefficient D s = linij^oo gz(t)/6t versus chain length N for different volume fractions <j>, as indicated in 
the figure (annealed SSA only). The large symbols correspond to the rescaled value ND B /Ab((j)) 2 where we use the acceptance 
rate A(<j>) ~ iV to characterize the local mobility. The small symbols (top) have been obtained by using the effective mobility 
m+(N,(/)), obtained from gi(t) (Fig. 0, as a measure of the local mobility instead of A. The dash-dotted line corresponds to 
the prediction in the dilute limit, = 2uo — 1 ft! 0.176. The data at <j) = 0.5 are compared to the effective power law with 
exponent (3 = —4/3 (dashed line) and the stretched exponential exp(— 0.8./V 1 / 3 ) (bold line). 
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FIG. 8: Rescaled times ti/ti (defined on page^J at <p — 0.5 (aSSA only). Similar behavior is also found for other densities. 
This confirms that all these times are proportional to one characteristic time scale — the "terminal time" Tt — containing the 
same information as the diffusion coefficients. The horizontal lines are comparisons with the Rouse model where t^/ti ~ 0.279, 
Ti/ti ~ 0.846 and ts/ti ~ 3.112. Here, we did not include r on d, as its scaling is different (see Fig.[5J. 
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FIG. 9: Unsealed times versus N. The open symbols present the results of the annealed SSA, the filled symbols of the quenched 
SSA. For the aSSA dynamics ti, r cn d and Td=3 (defined on pageQ are shown. In the dilute limit ((f) = 0.0078) we find n ~ TV 2 , 
as expected (small spheres, thin line). All other data are for (p = 0.5. The stretched exponentials are motivated by the ARH. 
They provide a good fit for both n w N 2 exp(0.8iV 1/3 ) (bold line) and r end « N 2/3v exp(0.8iV 1/3 ) (dash-dotted line) at high 
densities. For the qSSA dynamics the filled spheres, triangles and diamonds correspond to ri, t n / 2 (also defined on page|IJ 
and T+ q at <f> = 0.5. r+ q is the crossover time to free curvilinear diffusion for the qSSA (see Fig. 1101 for its definition in the 
aSSA case). As expected, the relaxation times for the qSSA increase much faster than for the aSSA. An effective power law 
t n/2 ~ iV 3 ' 5 is found in agreement with the corresponding diffusion coefficient D c (see Fig.EJ. The crossover to free curvilinear 
diffusion for the qSSA is roughly characterized by r+ q ~ iV 2 without noticeable curvature. 
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FIG. 10: Spatial MSD's for a dense system of chain length N = 512 at volume fraction (j> = 0.5 for the annealed SSA (pure 
slithering snake dynamics, i.e., uj = 0). The symbols are the same as in Fig. 2] While the behavior at large times is the same as 
in the previous figure, much slower dynamics and strong curvature are found at short times. This is a typical result for strongly 
overlapping chains (either large <f> or large N). According to Eq. Ill 21 we compare gi(t) = g2(t) « <?4 (t) = gs(t) and gs{t) with 
power laws characterized by an effective exponent a — 1/2. The time scale r+ characterizes the crossover between subdiffusive 
and free (linear slope) center of mass motion. Its scaling with N is discussed in Fig. 1141 




FIG. 11: Double logarithmic plot of spatial and curvilinear MSD's for the annealed SSA (open symbols) and the quenched 
SSA (filled symbols) for N = 128, <f> = 0.5 versus time. As suggested by Eq. {@}, we rescaled the curvilinear MSD ge(t) by 
Nl 2 /b 2 » 2N/3 to compare it with gs{t). For the aSSA spatial and curvilinear motion are coupled for all times: ga(t) ~ Ngs(t). 
For the qSSA the run is too short to reach the spatial free diffusion limit for N = 128. However, aSSA and qSSA dynamics 
are identical for small times t <C ta uc ~ r cn d (see Eq. C3J). Roughly logarithmic behavior is found for the qSSA dynamics, 
as indicated by the two fits gi(t) ~ log(t) and iV 173(4) ~ ge(t) ~ log(t) 2 . The time 7fl uc denotes the time where the density 
fluctuations of the annealed field become relevant, r coup the breakdown of the coupling of curvilinear and spatial displacements 
for the qSSA. 
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FIG. 12: Rescaled MSD of the curvilinear diffusion of the central monomer y — ge(t)/(2t). The MSD is normalized to yield 
the curvilinear diffusion coefficient D c as a plateau value. The annealed SSA data (open symbols) comprise all N for the 
density <j> = 0.5. For comparison, we also included the qSSA results for TV = 64, 128, 1024 (filled symbols). At short times 
(<T- k 50, all curves merge. For larger times the dynamics progressively slows down as a function of chain length. Only for 
times t S> t+(N) (indicated here for N = 1024) we obtain the expected free curvilinear diffusion. The intermediate decrease 
for the aSSA dynamics is very roughly fitted by the power law a — 1 » —1/2 (dashed line). The qSSA data for N — 1024 are 
well fitted by the logarithmic relation ge(t) 1 ^ 2 w log(t) (bold line). This suggests that for much longer chains even the aSSA 
dynamics might approach this logarithmic and chain length independent envelope. 
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FIG. 13: Relative scaling of the spatial and curvilinear displacements for the (a) aSSA (main figure) and (b) qSSA (inset) 
dynamics. The MSD's are defined on pageQ] (a) The spatial MSD's gi{t) and <?3(i) are plotted versus the curvilinear MSD 
ge (t) for various chain lengths, as indicated (4> — 0.5). As gs(t) is strictly linear with respect to ge,(t) for all t and N, we only 
included the data for N = 64 and N = 256 for the sake of clarity. As expected, we find gi(t) « b 2 (g 6 (t)/l 2 Y= 1/2 oc iV° for 
t C r t (dashed line labeled by v = 0.5). Hence, the anticipated relations between curvilinear and spatial diffusion are verified 
for the aSSA dynamics and the anomalous dynamics must be fully encapsulated in terms of the curvilinear displacement ge (t). 
(b) Rescaled spatial MSD 2Ngs(t)/3 versus g&(t). Only at short times and for densities <fi < 0.125 we find both quantities to be 
equal. The curvilinear motion decouples from the spatial motion for larger densities, where the snakes become localized within 
the correlation hole. 
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FIG. 14: Scaling attempts for ge(t) in order to determine the time scale at which the curvilinear dynamics (for the aSSA) 
becomes freely diffusive. The natural attempt r+ = ri (shifted for clarity) fails. In contrast, the scaling r+ = r cn d is successful, 
demonstrating that there is an additional characteristic time related to the density of chain ends. The power law with exponent 
a fa 1/2 for the envelope is indicated (dashed line). 




FIG. 15: Distribution of the first passage time tf p required to reach a curvilinear distance Nl/2. We consider aSSA systems at 
4> — 0.5 (main figure) and at chain length N = 256 (inset). Symbols are indicated in the figure. The data points of all systems 
(4>,N) collapse on the same master curve if rescaled by the mean first passage time Tf p = (tf p ). At large times the histograms 
decrease exponentially according to exp(— t/r c ), as expected. 
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FIG. 16: Normalized histograms of the curvilinear motion (aSSA data only) for N = 1024, <f> = 0.5 and for different time 
intervals At, as indicated in the figure. The histograms are rescaled in such a way that all data should collapse if they were 
self-similar for all times. We show that the histograms are approximately exponential within the time interval r_ w 50 -C t <C 
r+ = r cn d (AT = 1024) w 0.2 ■ 1Q 8 . The histograms become Gaussian again at much longer times. 
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FIG. 17: Non-gaussian parameter a.2 = | fyX ~ 1 for the curvilinear displacement s of the central monomer versus time for 
4> = 0.5 and chain lengths as indicated in the figure. We find that 02 becomes non-monotonous for r_ <(<t+. Unfortunately, 
the statistics is poor for large iV and some binning is necessary. Only the aSSA data are shown, the statistics is not sufficient 
for the qSSA. 
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FIG. 18: Spatial diffusion coefficient ND S versus N for different ui as indicated in the figure (aSSA only): u = corresponds 
to pure slithering snake dynamics, 1/lu = (large stars) to the pure local dynamics. The data for u<l are very similar to 
the pure slithering snake limit. For u ~ 8 Eq. @ is approximately satisfied. This may define a reasonable "starting point" for 
future SSA simulations corresponding to the classical URH reptation approach. For the local dynamics the power law prediction 
D B oc 1/N 2 (dashed line) from classical reptation theory is indicated. The data for uj = 512 approximate the local dynamics 
limit rather well for the chain lengths we have been able to simulate. 



